home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / nor.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  3KB  |  101 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. #define P10 242.66795523053175
  15. #define P11 21.979261618294152
  16. #define P12 6.9963834886191355
  17. #define P13 -.035609843701815385
  18. #define Q10 215.05887586986120
  19. #define Q11 91.164905404514901
  20. #define Q12 15.082797630407787
  21. #define Q13 1.0
  22.  
  23. #define P20 300.4592610201616005
  24. #define P21 451.9189537118729422
  25. #define P22 339.3208167343436870
  26. #define P23 152.9892850469404039
  27. #define P24 43.16222722205673530
  28. #define P25 7.211758250883093659
  29. #define P26 .5641955174789739711
  30. #define P27 -.0000001368648573827167067
  31. #define Q20 300.4592609569832933
  32. #define Q21 790.9509253278980272
  33. #define Q22 931.3540948506096211
  34. #define Q23 638.9802644656311665
  35. #define Q24 277.5854447439876434
  36. #define Q25 77.00015293522947295
  37. #define Q26 12.78272731962942351
  38. #define Q27 1.0
  39.  
  40. #define P30 -.00299610707703542174
  41. #define P31 -.0494730910623250734
  42. #define P32 -.226956593539686930
  43. #define P33 -.278661308609647788
  44. #define P34 -.0223192459734184686
  45. #define Q30 .0106209230528467918
  46. #define Q31 .191308926107829841
  47. #define Q32 1.05167510706793207
  48. #define Q33 1.98733201817135256
  49. #define Q34 1.0
  50.  
  51. #define SQRT2 1.414213562373095049
  52. #define SQRTPI 1.772453850905516027
  53.  
  54. void normbase(x, phi)
  55. double *x, *phi;
  56. {
  57.   int sn;
  58.   double R1, R2, y, y2, y3, y4, y5, y6, y7, erf, erfc, z, z2, z3, z4;
  59.  
  60.   y = *x / SQRT2;
  61.   if (y < 0) { 
  62.     y = -y;
  63.     sn = -1;
  64.   }
  65.   else sn = 1;
  66.   y2 = y * y;
  67.   y4 = y2 * y2;
  68.   y6 = y4 * y2;
  69.  
  70.   if(y < 0.46875) {
  71.     R1 = P10 + P11 * y2 + P12 * y4 + P13 * y6;
  72.     R2 = Q10 + Q11 * y2 + Q12 * y4 + Q13 * y6;
  73.     erf = y * R1 / R2;
  74.     if(sn == 1) *phi = 0.5 + 0.5*erf;
  75.     else *phi = 0.5 - 0.5*erf;
  76.   }
  77.   else if(y < 4.0) {
  78.     y3 = y2 * y;
  79.     y5 = y4 * y;
  80.     y7 = y6 * y;
  81.     R1 = P20 + P21 * y + P22 * y2 + P23 * y3 + P24 * y4 + P25 * y5 
  82.        + P26 * y6 + P27 * y7;
  83.     R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 + Q24 * y4 + Q25 * y5
  84.        + Q26 * y6 + Q27 * y7;
  85.     erfc = exp(-y2) * R1 / R2;
  86.     if(sn == 1) *phi = 1.0 - 0.5*erfc;
  87.     else *phi = 0.5*erfc;
  88.   }
  89.   else {
  90.     z = y4;
  91.     z2 = z * z;
  92.     z3 = z2 * z;
  93.     z4 = z2 * z2;
  94.     R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4;
  95.     R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4;
  96.     erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2));
  97.     if(sn == 1) *phi = 1.0 - 0.5*erfc;
  98.     else *phi = 0.5*erfc;
  99.   }
  100. }
  101.